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Abstract We prove sharp, computable error estimates for the propagation of er¬ 
rors in the numerical solution of ordinary differential equations. The new esti¬ 
mates extend previous estimates of the influence of data errors and discretisation 
errors with a new term accounting for the propagation of numerical round-off er¬ 
rors, showing that the accumulated round-off error is inversely proportional to the 
square root of the step size. As a consequence, the numeric precision eventually 
sets the limit for the pointwise computability of accurate solutions of any ODE. 
The theoretical results are supported by numerically computed solutions and error 
estimates for the Lorenz system and the van der Pol oscillator. 

Keywords High precision, high order, high accuracy, probabilistic error 
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1 Introduction 

We consider the numerical solution of general initial value problems for systems 
of ordinary differential equations (ODE), 

u{t) = t€{0,T], 

m ( 0 ) = uo, 

where the right-hand side / : x [0, T] —>• is assumed to be Lipschitz con- 

tinuous in u and continuous in t. Our objective is to analyse the error in an ap¬ 
proximate solution U : [0,T] computed by a single-step numerical method, 
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such as an explicit or implicit Runge-Kutta method. For the numerical results pre¬ 
sented at the end of this work, we have used a time-stepping method formulated 
as a Galerkin finite element method — which, for any particular choice of finite el¬ 
ement basis and quadrature, will correspond to a particular implicit Runge-Kutta 
method — but stress that the theoretical results are not particular to time-stepping 
based on Galerkin formulations but generic to all single-step methods. 

The propagation of local errors and accumulation of global errors in the numer¬ 
ical solution of ODE have been studied extensively in the literature, see e.g. 00 
010. These estimates are based on the formulation of an auxiliary dual problem: 
the linearised adjoint problem. From the solution of the dual problem, the accu¬ 
mulation rate of local errors may be computed, either as global stability factors 
or as local stability weights. These factors or weights, together with a measure of 
the local error, typically the residual R{t) = U — lead to a computable 

estimate of the global error. 

Standard estimates may include various sources contributing to the global 
error, such as discretisation errors, accounting for the use of finite time steps, 
quadrature errors, accounting for the approximation of the right-hand side / by a 
particular quadrature rule, and data errors, accounting for the approximation of 
the initial value uq. In this work, we extend these estimates by adding a new term 
accounting for the use of finite numeric precision in the computation of the numer¬ 
ical solution. This error is normally neglected, since it is typically much smaller 
than the contribution from the data or discretisation error. However, when the 
system o is very sensitive to perturbations, when the time interval [0,T] is very 
long, or when a solution is sought with very high accuracy, the effect of numerical 
round-off errors as a result of finite numeric precision can and will be the domi¬ 
nating error source, which ultimately limits the computability of a given problem. 


2 Main result 

We prove below that the error E in a computed numerical solution U approximat¬ 
ing the exact solution u of the ODE o is the sum of three contributions: 


E = E^) -t E(3 + Ec, 


where E^j is the data error, which is nonzero if U{0) i=- u(0); Eg is the discretisation 
error, which is nonzero as a result of a finite time step; and Eg is the computational 
error, which is nonzero as a result of finite numerical precision. Furthermore, we 
bound each of the three contributions as the product of a stability factor and a 
residual that measures the size of local contributions to the error. The size of the 
residuals may be estimated in terms of the size of the time step. We find that 

E - 5g(r)||t/(0) - u(0)|| + SG{T)Af + Sc{T)At-^l^, 

where At is the size of the time step, r the order of convergence of the numerical 
method, and Sg)(T), Sg(T), Sg(T) are stability factors which can be computed 
a posteriori. 
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3 Error analysis 

The error analysis is based on the solution of an auxiliary dual problem and follows 
the now well established techniques developed in ia,i6],i5] and [1] , with extensions 
to account for the accumulation of round-off errors. 

The dual (linearised adjoint) problem takes the form of an initial value problem 
for a system of linear ordinary differential equations: 

-i{t) = A^{t)z{t), te[0,T), 

z{T) = ZT. 


Here, A denotes the Jacobian matrix of the right-hand side / averaged over the 
approximate solution U and the exact solution u: 

= / ^{sU{t) + {1-s)u{t),t)ds. (3) 

By the chain rule, it follows that 

A{t){U{t) - u(t)) = j -h (1 - s)u{t),t){U{t) - u(t)) ds 

= j + (1 “ ds = - f{u{t),t). 

Based on the formulation of the dual problem we may now derive a (standard) 
error representation (Theorem [1]). It represents the error in an approximate so¬ 
lution U (computed by any numerical method) in terms of the residual R of the 
computed solution and the solution 2 of the dual problem ©■ The only assumption 
we make on the numerical solution U is that it is piecewise smooth on a partition 
of the interval [0,r] (or that it may be extended to such a function). At points 
where U is smooth, the residual is defined by 

R{t) = U{t)-f{U{t),t). 

Theorem 1 (Error representation) Letu : [0,r] — >■ be the exact solution of the 

initial value problem o, let z : [0,r] be the solution of the dual problem (f2|l . 

and let U : [0,T] —^ be any piecewise smooth approximation of u on a partition 
0 = to < ti <■■■< tM = T 0 / [0,r], that is, U\(^tm-i,t,n] ^ for 

m = 1,2,M {U is left-continuous). Then, the error U{T)—u{T) may be represented 
by 


M .T 

{zT, U{T) - u{T)) = (^(0), 1/(0) - m( 0)) + ^ (z(t™_i), [UU-A + / {z, R) dt 

m=l 


where R{t) = U{t) — f{U{t),t) is the residual of the approximate solution U and 
[U]m-1 = U{t+_A - U{tm-i) = U{t) - U{tm-i). 
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Proof By the definition of the dnal problem, we hnd that 

fT M tm 

(zT,e{T)) = (zT,e{T))- (i + e) dt = (zy, e(T)) - ^ / (i + A^z,e)dt, 

•^0 — l 


where e = U—u. Noting that {A^ z, e) = (z, Ae) and integrating by parts, we obtain 


M 


{zT,e{T)) = {z{0),em + 


{z{tm-i),[U]Tn-i) + {z,e-Ae)dt 


where J = C/(t+_ J -U{trn-i) denotes the jump of U 

at t = tm-i- By the construction of A, it follows that Ae = f{U, ■) — f{u, ■). Hence, 
e — Ae = U — f{U, ■) — u + f{u, ■) = U — f{U, ■) = R, which completes the proof. 

Remark 1 Theorem [T] holds for any piecewise smooth fnnction U : [0,T] —>■ 
in particular for any piecewise smooth extension of any approximate numerical 
solution obtained by any numerical method for ©■ 


We next investigate the contribution to the error in the computed numerical 
solution U from errors in initial data, numerical discretisation, and computation 
(round-off errors), E = -|- Eg + Ep, and derive sharp bounds for each term. 

To estimate the computational error, we introduce the discrete residual R 
defined as follows. For any p > 0, let be the Lagrange nodal basis 

for pP{[0,l]), the space of polynomials of degree < p on [0,1], on a partition 
0 < TO < n <•••< Tp < 1 of [0,1], that is, span{Afc}^^p = P^([0,1]) and 
Xi{Tj) = 5ij. Then, the discrete residual Rf^ is defined on each interval 
by 

RT = >^k{0)[U]m-l + [ ^ Xk{{t-tm-l)/Atn^)R{t)dt, fc = 0,l,...,p. (4) 

We also dehne the corresponding interpolation operator tt onto the space of piece- 
wise polynomial functions on the partition 0 = to < ti < • • • < = T by 

p 

{7Vv){t) = ^ v{tm-l + TkAtm) Xk{{t - tm-l)/Atm), t ^ [tm-lAm]- 
k=0 

We may now prove the following a posteriori error estimate. 

Theorem 2 (Error estimate) Let u : [0,r] —>■ he the exact solution of the 

initial value problem let z : [0,r] —>■ be the solution of the dual problem ([2)1 . 

and let U : [0, T] —>■ be any piecewise smooth approximation of u on a partition 

0 = to < ti < ■■■ < tM = T of [0,T], that is, ^ C°°((tm-i,tmj) for 

m = 1,2, ...,M {U is left-continuous). Then, for any p > 0, the following error 
estimate holds: 

{zT,UiT)-u{T)) = ED + EG + Ec, (5) 

where 

lEiol <S,3||t/(0)-u(0)||, 

|Eg| < SGCpmax{z\tP+i(||[[/]||/4it+ ||7?||)} , 

|Ecl < SnCL max max 
' Po<k<p[o.T]" " 
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Here, Cp and Cp are constants depending only onp. The stability factors So, Sq, and 
Sc are defined by 


Sd = l|2 


I 

L 


Sg= I 


Sc = / l|7r-*l|df. 

Jo 

Proof Starting from the error representation of Theorem [1] we add and subtract 
the degree p left-continuous piecewise polynomial interpolant ttz defined above to 
obtain 


(zT,e{T)) = (z(0),e(0)) 


M 


+ E 

m=l 

M 

+ E 


m=l . 

= E£) -I- E(3 + Ep. 


{z{tm-l) - r:z{tl^_^),[U]m-l) + [ {z-'kz,R) 

Jtm-l 

(7r2(t+_i), [(7]m-i) + / {T:z,R)dt 


dt 


We first note that the data error E^j is bounded by ||2(0)|| ||e(0)|| = So ||e(0)||. By 
an interpolation estimate, we may estimate the discretisation error Eg by 


M 


z{tm-l) -T^z{t:^_^)\\\\[U]m-l\\+ [ \\z-nz\\\\R\\dt 

J tjfl _1 

M O. 


Eg < E 

m=l _ 

Cp max {AtP+\\\[U]\\/At+\\R\\)}j2 r 

^ ’ ‘ m,= 1 dtm-1 


where ^ dt = Jq dt = Sc and Cp is an interpolation 

constant. Finally, to estimate the computational error, we expand ttz in the nodal 
basis to obtain 


Eg = 


M P / \ 

EE z{tm-l+'rkZ^tm),\k{d)[U]m-l+ / \k{{t - tm-l) / Atm)R{t) dt\ 
m.= lk=o\ Jim -1 / 


m=i k= 

M p M p 

= ^ '^{z{tm-l + rf,At„i), RT) = E ^trn'^{z{tm-l + TkAtm), Atfi,^ R'^) 

m=l k=0 m=l k=0 

M p 

< E E WAtfn^R'^W 

m=l k=0 

M P 

< max max||Zit“^.Rj,|| Atmy^ \\z{tm-i + rkAtm)\\ 

0<k<p[0,T] y 

- - ^ ^ m=l k=0 

M rtn 


< CL max max IIZ\t ^Ri-W 

- ^0<k<p[0,T]" 


E/" 

m=l -'im-i 


\ttz\\ dt, 
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where J2m=i It"' i ~ lo = '^C' and Cp is a constant depending only 

on p. This completes the proof. 

Remark 2 Theorem[2]estimates the size of {zp, U{T)—u{T)) for any given vector zp- 
We may thus estimate any bounded linear functional of the error at the hnal time 
by choosing zp as the corresponding Riesz representer. In particular, we may 
estimate the error in any component Ui {T) of the solution by setting zp to the ith 
unit vector for i = 1,2,..., A'". 


Theorem[2]extends standard a posteriori error estimates for systems of ordinary 
differential equations in two ways. First, it does not make any assumption on the 
underlying numerical method. Second, it includes the effect of numerical round-off 
errors. A similar estimate can be found in m but only for the simplest case of 
the piecewise linear cG(l) method (Crank-Nicolson). 

We now investigate the propagation of numerical round-off errors in more de¬ 
tail. As in the proof of Theorem [2l Eq denotes the computational error dehned 
by 


M 




(7r2{t+-l)> [Ujm-l) + 


L 


{■RZ, R) dt 


( 6 ) 


Theorem[2]bounds the computational error in terms of the discrete residual defined 
in &■ The discrete residual tests the continuous residual R = U — f oi dlD against 
polynomials of degree p. In particular, it tests how well the numerical method 
satisfies the relation 


U{tm) 



/{[/,■) dt. 


( 7 ) 


With a machine precision of size Cmachi our best hope is that the numerical method 
satisfies © to within a tolerance of size Cmach for each component of the vector U. 
It follows by the Cauchy-Schwarz inequality that 

max HR™ 11 < Cmach^iV. 
k,m 


We thus have the following corollary. 


Corollary 1 The computational error Ec of Theorem\^is bounded by 


|Ecl < Sc 


^mach^ R 
^ minjQ At 


This indicates that the computational error scales like At~^\ the smaller the time 
step, the larger the computational error. At first, this seems non-intuitive, but it is 
a simple consequence of the fact that a smaller time step leads to a larger number 
of time steps and thus a larger number of round-off errors. 

The estimate of Corollary [1] is overly pessimistic. It is based on the assumption 
that round-off errors accumulate without cancellation. In practice, the round-off 
error is sometimes positive and sometimes negative. As a simple model, we make 
the assumption that the round-off error is a random variable which takes the value 
+fmach or “^mach with equal probabilities. 



“bCmachi 

^mach? 


P = 0.5, 
p = 0.5. 


(8) 
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for all m, k, i. In reality, round-off errors are not uncorrelated random variables, 
but the simple model m may still give useful results. For a discussion on the 
applicability of random models to the propagation of round-off errors, see m 
(Section 2.8) and [TO] . 

Under the assumption m, we find that the expected size of the computational 
error scales like . As we shall see in the next section, this is also confirmed 

by numerical experiments. A similar result is obtained in a series of papers by 
Li et. al |191I20| . In [T^, it is first noted that there exists an optimal time step-, that 
is, a time step for which discretisation errors and round-off errors balance. In |20| . 
it is then found that the round-off error is inversely proportional to the square 
root of the time step. These results are confirmed by the following theorem. 


Theorem 3 Assume that the round-off error is a random variable of size 

with equal probabilities. Then the root-mean squared expected computational error Ec 

of Theorem\^ is bounded by 

(E[E^])1/2 < Sc, 

'' min[Q (p] y/ At 


( cT Ii2 1 \ / 

where Sc^ ~ [Jo ll^'^ll ^ constant depending only onp. 

Proof As in the proof of Theorem [2l we obtain 

M p M p N 

Ec=EE {z(tm-l + TkAtm), RT) = EEE Zi{tm—1 T 'Ck^l^rn)[Rk )*! 

m=lfe=0 m=lk=0i=l 

where by assumption {R^)i = ejnach^mfci ^mki — with probability 0.5 and 
0.5, respectively. It follows that 

M p N 

E E E [tm— 1 + ri;Atm)zj(tn — l + TiAtn) Cmach^ mki^nlj 

m,n=l k,l=0 i,j=l 

^ ^ (fm — l ~\~ ri^Atm) ^mach^mfci 

{m,k,i) = {n,l,j) 

-h E m —1 + TfeAtm)2j(tn-l -\-TiAtn) CmaLch^mki^nlj • 

We now note that = 1. Furthermore, yijkimn = ^mki^nij is a random variable 
which takes the values -|-1 and —1 with equal probabilities. We thus find that 

M p N 

El^'c] = ^mach EEE zfitm-1 + TkAtm) + 0 

m=l k=Q i=l 
M p 


^mach 


E E 


m=l k=0 
M 


< E E 1 + < Sl^ C'p 

mm\n At ^mm[o,T]^^ 


” m=l fc=0 
vl/2 


/ T 2 

where Sc, = (/g ||7r2|pdt j . This completes the proof. 
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Remark 3 By additional assumptions on the smoothness of the dual solution z, one 
may relate the error to the expected value of the distance from the starting point for 
a random walk which is y/n for n steps (for n large), and prove a similar estimate for 
the expected absolute value of the computational error, i?[|Ec|] ~ ScCmach/V^St, 
where Sc = Jq ||7r-z|| dt. 

Remark 4 By Cauchy-Schwarz, the stability factor Sc of Theorem [5] is bounded 
by y/TSc^- 

Remark 5 In m, the effect of numerical round-off error accumulation and its rela¬ 
tion to Brownian motion (Brouwer’s law) are discussed in the context of symplectic 
methods for Hamiltonian systems. It should be noted that although the assump¬ 
tions of Theorem |3] are similar to those in m, namely that the process of error 
accumulation for round-off errors is random rather than systematic, the point un¬ 
der discussion in the present work is different: the effect of time step size rather 
than the effect of the interval length. 


We conclude this section by discussing how the above error estimates apply to 
the particular methods used in this work. The estimate of Theorem [2] is valid for 
any numerical method but is of particular interest as an a posteriori error estimate 
for the hnite element methods cG{q) and dG(q) (see |121ll3l[m il]l. 

The continuous and discontinuous Galerkin methods cG{q) and dG(q) are for¬ 
mulated by requiring that the residual R = U — f(U, •) be orthogonal to a suitable 
space of test functions. By making a piecewise polynomial ansatz, the solution 
may be computed on a sequence of intervals partitioning the computational do¬ 
main [0, T] by solving a system of equations for the degrees of freedom on each 
consecutive interval. For a particular choice of numerical quadrature and degree q, 
the cG( 5 ) and dG((/) methods both reduce to standard implicit Runge-Kutta 
methods. 

In the case of the cG((ir) method, the numerical solution [/ is a continuous 
piecewise polynomial of degree q that on each interval satishes 



V Rdt = Q 


( 9 ) 


for all V € It follows that the discrete residual Q is zero if p < 

(jr—1. However, this is only true in exact arithmetic. In practice, the discrete residual 
is nonzero and measures how well we solve the cG{q) equations (|9l), including 
round-off errors and errors from numerical quadrature^ For the cG{q) method, 
we further expect the residual to converge as At‘^. Thus, choosing p = q — I in 
Theorem [2l one may expect the error for the cG(g) method to scale as 

F = E£) -|- Eq + E(7 ~ S{T) ^Emach + + At ■ (10) 

Here, S{T) denotes a generic stability factor. As in Theorem [2] each term con¬ 
tributing to the total error is in reality multiplied by a particular stability factor. 
In practice, however, the growth rates of the different stability factors are similar 
and related by a constant factor. 


^ To account for additional quadrature errors present if the integral of IS} is approximated 
by quadrature, one may add and subtract an interpolant tt/ of the right-hand side / in the 
proof of Theorem [ 2 ] to obtain an additional term Eg = 5 q max[o,T] IK/ ~ /II where Sq = 

lo ~ ^C- 
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4 Numerical results 

In this section, we present numerical results in support of Theorem [2] and Theo- 
rem[3l The examples are the well-known Lorenz system and Van der Pol oscillator. 
Both examples illustrate the competing convergence rates for discretisation errors, 
decreasing rapidly for smaller time steps, and computational errors (round-off er¬ 
ror), increasing for smaller time steps. 

The numerical results were obtained using the authors’ software package Tan¬ 
ganyika |23| which implements the methods described in [21] using high precision 
numerics provided by GMP [9]. A complete code for reproducing all results in this 
paper is available at m- For details on the implementation, see [16] . 


4.1 The Lorenz system 

We first consider the well-known Lorenz system [23], a simple system of three 
ordinary differential equations exhibiting rapid amplification of numerical errors: 

' x = a(y-x), 

< y = rx — y — xz, (11) 

z = xy — bz, 

where ct = 10, fo = 8/3, and r = 28. We take u(0) = (1,0,0). 

The Lorenz system is deterministically chaotic. In the context of a posteriori 
error analysis of numerical methods for the solution of ODE initial value problems, 
as in the present work, this means that solutions may, in principle, be computed 
over arbitrarily long time intervals, but to a rapidly increasing cost as function of 
the final time T. 

4-1.1 Computability and growth of stability factors 

In [8], computability was demonstrated and quantified for the Lorenz on time 
intervals of moderate length (T = 30) on a standard desktop computer. This result 
was further extended to time T = 48 in [52], using high order (||e(T)|| ~ At^^) finite 
element methods. Solutions over longer time intervals have been computed based 
on shadowing (the existence of a nearby exact solution), see [3], but for unknown 
initial data. Related work on high-precision numerical methods applied to the 
Lorenz system include [^ and m- 

In [TJ, the authors study the computability of the Lorenz system in detail 

on the time interval [0,1000]. Computability is here defined as the maximal final 

time T = r(emach) such that a solution may be computed with a given machine 

precision tmach- The computability may be estimated by examining the growth rate 

of the stability factors appearing in the error estimate of Theorem[2] By numerical 

solution of the dual problem, it was found in m that the stability factors grow 

exponentially as S{T) ~ ^qO-SSST ^ 10° "*'^; see Figure [T] By examining in detail 

the terms contributing to the error estimate m, one finds that an optimal step 
1 

size is given by At ~ c^a'ch that the computability of the Lorenz system is 
given by 

T(Cinach) ^ 
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T T 

Fig. 1 Growth of the stability factor Sc (left) for the Lorenz system on the time inter¬ 
val [0,1000] and a detailed plot on the time interval [0, 50] (right). 


where rimach = ~ logic ^mach Is the number of significant digits. Based on this 
estimate, one may conclude that with 16-digit precision, the Lorenz system is 
computable on [0,40], while using 400 digits, the Lorenz system is computable on 
[ 0 , 1000 ]. 

In Figure[2l we plot the solution of the Lorenz system on the interval [0,1000]. 
The solution was computed with cG(lOO), which is a method of order 2q = 200, 
a time step of size At = 0.0037, 420-digit precision arithmetic^, and a tolerance 
for the discrete residual of size tmach ~ 2.26 • 10“^^^. The very rapid (exponen¬ 
tial) accumulation of numerical errors makes the Lorenz “fingerprint” displayed in 
Figure [2] useful as a reference for verification of solutions of the Lorenz system. If 
a solution is only slightly wrong, the error is quickly magnified so that the error 
becomes visible by a direct inspection of a plot of the solution. 


4-1.2 Order of convergence and optimal step size 

We next investigate how the accumulated error at final time depends on the size 
of the time step At. According to (HOD, we expect the error to scale like At^‘‘ + 
Thus, for a gradually decreasing step size, we expect the error to 
decrease at a rate of At'^'^. However, as the time step becomes smaller the second 
term At~^^^ will grow and, for small enough At, be the dominating contribution 
to the error. This picture is confirmed by the results presented in Figure[3]for two 
numerical methods, the 2nd order cG(l) and the 10th order cG(5) method. Of 
particular interest in this figure is the very short range in which the 10th order 
convergence of the cG(5) method is recovered; with only 16 digits of precision, the 
dominating contribution to the total error is the accumulated round-off error. We 
also note that for both methods, one may find an optimal size of the time step At 
for which both contributions to the total error are balanced. 

In Figure m results are presented for an investigation of the influence on both 
the step size At and the polynomial degree g in the cG(g) method. As expected, 
the minimal error is obtained when both the polynomial degree g and the step 

^ The requested precision from GMP was 420 digits. The actual precision is somewhat higher 
depending on the number of significant bits chosen by GMP. 
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,0 I-^^^^-1 

900 920 940 960 980 1000 


t 

Fig. 2 Accurate reference solution for the three components of the Lorenz system on the 
interval [0,1000] with the x and y components plotted in blue and green respectively (and 
almost overlaid) and the z component in red. 
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Fig. 3 Error at time T = 30 for the cG(l) solution (left) and at time T = 40 for the cG(5) 
solution (right) of the Lorenz system. The slopes of the green lines are —0.35 ~ —1/2 and 
1.95 ~ 2 for the cG(l) method. For the cG(5) method, the slopes are —0.49 ~ —1/2 and 
10.00 Ri 10. 


size are maximal. Maximising the step size minimizes the influence of numerical 
round-off errors (the term and as a consequence the polynomial degree q 

must be large in order to suppress the discretisation error (the term 


4.2 The Van der Pol oscillator 

We next consider the Van der Pol oscillator, given by the second order ODE 

it = /t(l — M^)it — u. 


Rewritten as a system of first order equations, it reads 


ill = U2, 

U2 = At(l — u\)u2 — tti 

We compute solutions on [0,2/i] for /r = 10^ and u(0) = (2,0). This configuration 
is used as a test problem for ODE solvers in [25]. Eor large values of the parameter 
/i, the solution quickly approaches a limit cycle. 

As for the Lorenz system, the stability factors(s) grow very rapidly (exponen¬ 
tially), as indicated in Eigures |5| and |6| However, the rapid growth is localised 
in time close to T « 807 • n for n G N. For times before or after these points of 
instability, the stability factor is of moderate size. This means that solutions are 
difficult to compute only at points near the points of instability; that is, a solution 
may be easily computed at time t = 1000 but not at time t = 807. 

This rapid growth of stability factors is reflected in the growth of the error 
for numerical solutions as shown in Figures |7| and |8| Examining these plots in 
more detail, we notice that in accordance with the error estimate of Theorem [2] 
and Equation (HQi). For the numerical solutions studied in Figures |7| and jS] the 
discretisation error dominates for low order methods. As the polynomial degree 
q is increased, the error decreases until the point when the computational error 
starts to dominate. We notice that the baseline error is of size E ~ for 


(13) 
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Fig. 4 Top: The accumulated total error at final time T = 40 for numerical solutions of the 
Lorenz system with different step size At and polynomial degree q using the cG(g) method. 
Due to the random nature of the round-off errors, the data has been smoothed. Lower left: 
Contour lines of the smoothed data. Lower right: The raw data included for completeness. 


the highest order methods when the stability factor is of size S ~ 10^, and the 
error spikes a,t E ^ 10“^ at times when the stability factor takes on large values 
S ^ 10^^. This is in good agreement with the error estimate: E ^ S ■ 
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Fig. 5 Growth of the computational stability factor Sc{'^) for fh® Van der Pol oscillator (USD . 




Fig. 6 Detail of growth of the computational stability factor Sc{T) for the Van der Pol 
oscillator ITSt . 


5 Conclusions 

We have proved error estimates accounting for data, discretisation and compu¬ 
tational (round-off) errors in the numerical solution of initial value problems for 
ordinary differential equations. These error estimates quantify the accumulation 
rates for numerical round-off error as inversely proportional to the square root of 
the step size, and proportional to a specific computable stability factor. The effect 
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Fig. 7 Growth of error for solutions of the Van der Pol oscillator II13II computed with time 
step At = 10-3. 




Fig. 8 Detail of growth of error for solutions of the Van der Pol oscillator II13II computed with 
time step At = IQ—®. 


of round-off errors is mostly pronounced for large values of the stability factor, 
which includes both chaotic dynamical systems as well as long-time integration of 
systems which exhibit only a moderate growth of the stability factor. 
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